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We apply optimum experimental design (OED) to organic semiconductors modeled by the extended Gaussian 
disorder model (EGDM) which was developed by Pasveer et al. 1 We present an extended Gummel method to 
decouple the corresponding system of equations and use automatic differentiation to get derivatives with the 
required accuracy for OED. We show in two examples, whose parameters are taken from Pasveer et ali and 
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I. INTRODUCTION 

Simulation of organic semiconductor devices, e.g. or- 
ganic light emitting diodes (OLED), organic solar cells, 
etc., has gained great interest in the past decade. Ac- 
curate models often lack in precise parameters and their 
identification is time-consuming or expensive. One ma- 
jor problem is the uncertainty in the measurement data 
which leads to uncertain parameters. Carrying out more 
experiments would minimize this uncertainty in an ex- 
pensive way. An alternative is to use OED in order 
to minimize the parameter uncertainty by planning new 
(optimal) experiments. New measurement data are re- 
ceived for which the parameter estimation yields param- 
eters with minimal confidence intervals. We apply the 
concept of OED to the EGDM, a special model for the 
mobility of electron transport in organic polymeric ma- 
terial. The chapters are arranged as follows: We give 
a brief overview of the model equations in Sec. [XT] and 
point out what the relevant quantities are. In Sec. IIII| 
we explain the methodology of the optimum experimen- 
tal design problem in more detail. After describing the 
equation solver methods we used, Sec. IIV[ numerical re- 
sults of the optimization are presented in Sec. [Vj 



II. EGDM EQUATIONS 

A basic description of charge transport in semiconduct- 
ing materials in the steady-state case is given by the cou- 
pled van Roosbroeck system 3 consisting of the continu- 
ity equation, also called drift-diffusion equation, and the 
Poisson equation. Given a domain Q := (0,L) C R, the 
state variables, i.e. the space dependent functions, are 
the electric charge density n in [m~ 3 ] and the electric 



potential <j) in [eV] which are scalar real valued functions 
defined on Q. We assume that n and <\> are twice dif- 
ferentiate on Q. Pasveer et al. 1 proposed the EGDM 
to conjugated semi-conducting polymers, where the dif- 
fusion and the mobility depend on the state variables 
and hence on the space variable x. Furthermore they in- 
troduced another state Ep, called Quasi- Fermi-Energy, 
and a corresponding equation which couples Ep and n 
at every space point. We omit the ^-dependence in the 
following equations, only n, <ft and Ep are space depen- 
dent. The model equations are: 



= d x J(n,(/>, E F ), 
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In these equations J is the electric current density in 
[-4r]> ks is the Boltzmann constant in [^jf], £ the permit- 
tivity in [y^], e the elementary charge in [As], a := 



2 



and 5(a) := 2- 



clog 4 



In the anorganic case, 
the gi-factors, z = 1,2,3, would be constant. Their or- 
ganic model is yield by comparison with the solution of 
the master equation. 1 On the boundary dft = {0,L} the 
following conditions are imposed: 
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where V is the voltage in [V] and y?i ? 2 are given energy 
barriers in [eV]. On the entrance, x = 0, the energy 
barrier y?i is lowered according to the theory of Emtage 
and O'Dwyer— and Scott and Malliaras£ For our compu- 
tations we take the dimensionless form of the equations 
proposed by Bonham and Jarvis. 6 For the later use we 
define 



p := (fio,a,N t ) G 
q:= (L,T)gR 2 , 



(11) 



where p are parameters of unknown numerical value given 
by nature. They have to be identified by comparing a 
model response to experimental data, /xq is the zero tem- 
perature mobility in [yj], a is the width of the Gaussian 
distribution of the density of states in [eV] and N t is the 
site density in [ra -3 ]. We assemble quantities which are 
adjustable by an experimenter, the device length L in [m] 
and the temperature T in [K], in the vector qi 



III. SOLUTION METHODS FOR THE EGDM 
EQUATIONS 

There are two established ways for solving the system 

1. Apply Newton's method to the fully coupled sys- 
tem of equations 

2. Use Gummel's method, 9 i.e. a fixed point iteration 
which decouples the three EGDM equations 

With Newton's method, one can achieve quadratic con- 
vergence. However, finding good starting values is not 
simple. In the work of Knapp et al., 8 a strategy moti- 
vated by physical considerations is described. Another 
aspect is that the Jacobian, i.e. derivatives of the func- 
tions w.r.t. the states, is required. To compute the Jaco- 
bian, the main two options are to use difference formulas 
or compute the exact derivatives. The latter can be done 



TABLE I. Algorithm for the fixed point iteration of Gum- 
mel expanded by the EGDM-Quasi-Fermi-Equation and a 
derivative- free linearization. 

Let u° := (n°,(jP,E%) be given and choose ||Au°|| > TOL 
with given error tolerance TOL. 
Set i = 0. 

while ||A^|| 2 > TOL 
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for E 7 ^ 1 with Newton's method started with E % F . 
Solve 

-dU i+1 = \n x 
with corresponding boundary conditions for 
With n\ and E^ 1 solve 

O = d x {g 1 (n i )g 2 (0 i+1 ) 

(n i+1 d x <P +1 - k B Tg 3 (n\ d x n i+1 

with corresponding boundary conditions for n z+1 . 
Set u i+1 := (n l+1 ,^ +1 ,^ +1 ) and Au t+1 := u i+1 - 



by hand or by using automatic differentiation (AD). Ei- 
ther way, the additional effort for computing the deriva- 
tive in n directions is at least 2n times the effort of each 
function evaluation in every Newton step. Another pos- 
sibility is to solve ([I])-© with Gummel's method. TA- 
BLE [J shows the modified algorithm consisting of the 
classical system of equations (pQ) and ([2]) and the Quasi- 
Fermi energy Ep defining equation ((3]). 10 We solve the 
equations sequentially and insert the interim results into 
the next equations. With the special linearization, used 
in the third equation of TABLE HI we do not need any 
derivatives and the effort of function evaluations is lim- 
ited. We discretize the infinite dimensional problem (pQ)- 
(|2j) to a finite one. The domain O is divided in N subin- 
tervals 1{ := [xi, Xi + h], X{ G {xq = 0, . . . , xn = L} with 
a constant mesh size h := ^. Finite differences are ap- 
plied to the spatial derivatives, i.e. with respect to x. The 
so-called Scharfetter-Gummel scheme^ forces the func- 
tion J, defined in Eq.(|5]), to be constant on each interval 
Ii, denoted by and provides an upwind stabiliza- 

tion, so that computation on coarse meshes is possible. 
On the interval L, the scheme looks like 



■ n i+i exp 
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The terms g^ j = 1,2,3 stand for average values of the 
non-constant functions gj, j = 1,2,3. It is important 
that the averages are taken in an upwind conform way to 
prevent numerical oscillations. 
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IV. OPTIMUM EXPERIMENTAL DESIGN FOR MODEL 
VALIDATION 

In this part, we follow the approaches of Lohmann 12 
and Korkel et al^ With different choices of the controls 
q and the voltage V, defined in Sec. Ull we set up multiple 
experiments in which the current density J ([Sj) is mea- 
sured. Let M be the number of measurements we yield. 
In a parameter estimation, the parameters are identified 
by fitting a model response, here J, to experimental data, 
i.e. measurements. We assume the measurement error to 
be normally distributed with mean zero and covariance 
matrix S 2 = diag(a 2 ,z = 1,...,M) G R MxM . With 
the same experimental settings, i.e. equal controls g, a 
fit from a different realization of the measurement errors 
may result in very different parameter values. The co- 
variance matrix of the parameters allows to analyze the 
quality of a parameter estimation. The assumed model 
for the standard deviations of the measurement errors 

IS 3^ 



with Si = (Axj)j and £2 = (Ayk)k- We define the 



Vi = 0.1 • Ji +0.1 



A 



where Ji G M is the function value of J corresponding to 
the z-th measurement. For further notation, we assemble 
the values Ji in the vector J G M M . If the confidence 
region of the parameters is approximated by assuming a 
linear propagation of the measurement errors, it can be 
parameterized by the covariance matrix defined by 



Cov p :=E[(p-E[p])(p-E[p}) T ] G 



t>N v xN v 



where N p is the number of parameters. From now on, we 
denote by n, <\> and Ep the discrete counterparts of the 
state variables which are (N — l)-dimensional vectors, 
without boundary values. They assemble the function 
values at the mesh points o^, i = 1, . . . , N — 1, cf. Sec JIIII 
We abbreviate the discrete solution of the system ([I])-(j3]) 
dependent on parameters p and controls q by 



u(p,q) := (n(p,q),</)(p,q),E F (p,q)) G 



cf. TABLE ID We used the notation N u := 3(7V - 1) for 
the overall state dimension. In the following we denote 
the derivative of a function / w.r.t. x in the direction Ax 
by 



d x f{Ax} := lim 



f(x + eAx) 



and write accordingly the second derivatives of / w.r.t. 
x and y in the directions Ax and Ay as dy X f{Ax, Ay}. 
We combine several directions Axi into a so-called seed 
matrix S — (Axi, • • • , Ax n ) and define d x f{S} by 

(d x f{S}) tJ = dxf^xj} 

and accordingly d% x f{Si, S2} by 

(d 2 yx f{S 1 ,S 2 }) ijk = dlJiiAxj^Vk} 



M x N p Matrix 



Jac(u(p,q),p,q) := -S 1 d p J(u(p,q),p,q){I p } 



-s- 1 



d u j{d p u{i p }} + d p j{i p } 



with the derivative of J w.r.t. to the parameters p with 
the TVp-dimensional identity matrix I p as seed matrix. 
Computing d u J{d p u{I p }} is much less expensive than 
computing the matrix product d u J{I u } ■ d p u{I p } with 
the identity matrix I u G M. Nu x Nu . The covariance matrix 
in the unconstrained case can be computed by 

Cov p = ( J&c(u(p,q),p,q) T Jac(u(p,q),p,q)) \ (12) 

For given probability a G [0, 1] the linearized (100 • a)%- 
confidence region is described by 

G(a,p) = {ve R Np : v = p + 6p, Sp T Gov' 1 Sp < 7(a) 2 }, 

(13) 

with the (1 — a)-quantile of the x 2_ distribution 7(a) 2 . 
As an approximation of the confidence intervals of the 
parameters we use 



\pi-0i,Pi + 9i] with 



<y(a)^(Cov p )ii, i = l,...,iV p . 



(14) 



We point out that for computing Cov p (and the confi- 
dence intervals) no measurement data is necessary. We 
consider the reduced nonlinear optimization problem 



1 



min — — trace ( Cov p 



(15) 



which is called "optimum experimental design problem" . 
Efficient optimization algorithms require derivative infor- 
mation. For higher order difference formulas, the number 
of function evaluations increases which can be very costly 
like in our case. Additionally, low order schemes reduce 
the number of exact digits, in particular in combination 
with derivatives of higher order. To avoid such approxi- 
mation errors, we use automatic differentiation where the 
derivatives are evaluated up to machine precision. For 
derivative based optimization, we need the derivative of 
the gradient of Jac with the seed matrix I q G M NqXNq , 
where N q is the number of the controls. N p and N q are 
relatively small, so we use the forward mode, where we 
have to compute 



d q Jac(ii(p, q),p, q){I q } 

~d q (d u J{d p u{I p }} + d p J{I p }){I q } 



%J{d P u{I p }, d q u{I q }} + d 2 qu J{d p u{I p }, I q } 
+ d u J{d 2 qp u{I pj I q } } + d 2 up J{I P , d q u{I q }} 
+ d 2 qp J{I p J q } 
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see e.g. Griewank^ We need derivatives of u w.r.t. the 
parameters controls q and mixed second derivatives 
w.r.t. p and q. One way would be to differentiate the fixed 
point iteration, TABLE HI and iterate over the derivatives 
as well. 15 However, we are only interested in the deriva- 
tives of the solution of Gummel's method u{p,q). In the 
solution it holds 

F(u(p,q),p,q) = Vp, q, (16) 

with F the discrete counterpart of the system of equa- 
tions (pQ)-(|3]). To compute the required derivatives 
d p u{I p }, d q u{I q } and d^ p u{l pi I q }, we use the sensitivity 
method, see e.g. Hinze et alJ£ Using adjoint method is 
not recommended here, because the number of measure- 
ments M is much higher than the number of parameters 
N p . Differentiating (fTSj) leads to 

d p F(u(p, q),p, q){I p } = 0, 
d q F(u(p, q),p, q){I q } = 0, 
d 2 q pF(u(p,q),p,q){Ip,I q } = 0. 

Hence the required derivatives are given in form of 

d pU {i p } = -(a„j , {j u })" 1 a p F{/ p } 

d qU {I g } = -{duFil^dqFilg} 

d 2 qp u{I p , I q } = - (duFih})- 1 (d 2 qp F{I p , I q } (17) 
+ 8 2 up F {I p , d q u{I q }} + 3 2 qu F {8 p u{I p }, JJ 

+ d 2 u F{d pU {i p },d pU {i q }}). 

We can save the decomposition of the matrix d u F{I u } 
and use it for all three equations to reduce complexity. 

V. NUMERICAL RESULTS 

Experimenters proceed in the following way^ Devices 
of different lengths L are produced by evaporating. At 
different temperatures T, they apply voltage series V, 
a vector consisting of various voltages V, and measure 
the corresponding current values, combined in the vector 
J . One choice of L, T and voltage series V we call one 
experiment subsequently. We will show in two examples 
how the confidence regions ([T3]) can be reduced by OED. 
In each example we consider three different lengths 

£> '= (Lij L 2 , Ls) 

and three different temperatures 

T:= (T l5 T 2 ,T 3 ). 

The combination of each element of the vector C with 
each element of the vector T lead to nine experiments. 
We assemble all control variables, which are the opti- 
mization variables for the OED problem, in a vector 

Q:= (L 1 ,T 1 ,L 1 ,T 2 ,L 1 ,T3,L 2 ,T 2 ...,L 3 ,T 3 ) G R 18 . 



The result of the OED should be three lengths and three 
temperatures combined in nine experiments. To enforce 
that, we apply additional constraints to Q for the lengths 

Si = Qs = e 5 , 
er = e 9 = Sii, 

2l3 = 2l5 = Ql7? 

and for the temperatures 

22 = 28 = 2l4, 
Q4 = 2l0 = 2l6; 

2 6 = 2i2 = 2i 8 . 

We solve problem ([T5]) with the software package 
VPLAN 18 where the optimization problem is solved with 
an inexact SQP method provided by SNOPT 7.2-9 (Jun 
2008). In each iteration, the objective, i.e. the covari- 
ance matrix, is computed by solving the underlying sys- 
tem of equations. We implemented the extended Gum- 
mel method TABLE U in VPLAN to solve this system. 
For the optimization, the derivatives of the model func- 
tions are calculated with the AD tool ADIFOR 2.O. 20 
We implemented a routine in VPLAN, which solves the 
system of equations (fT7|) . In a first example we take con- 
figurations p, V, C and T from Pasveer et ali who used 

p = (1.15 • 10" 5 , 0.14, 2.44 • 10 26 ) , 
V = (0.1,..., 0.9, 1,2,..., 10), 
C = (275,350,475), 
T = (235,270,305). 

Here and in the following example we set cpi = </? 2 = 
for the energy barriers in [eV] and e = 2.66 • 10 -11 for 
the permittivity in [y^]- We choose boundaries to the 
lengths and temperatures according to the practicability 
of experiments and the validity of the model: 

Li e [50, 500] and T { e [200, 350] i = 1, 2, 3. (19) 

The optimization results in the vectors 

C* = (50,339.1,471.6), 
T* = (277.2,281.8,350). 

Where the lengths do not change much, a tendency to 
higher temperatures is observed. The algorithm gives 
back the exact numbers of the optimum. Dependent on 
the equipment of the laboratory these values can not be 
realized in practice. Nevertheless they can be taken as a 
guideline. All the values in the neighborhood of C* and 
T* are leading to similar results. In Fig. [TJ we visu- 
alize the three dimensional ellipsoid and the projections 
corresponding to the set Gl(0.95, 0) before and after the 
optimization. TABLE HI1 shows a comparison of the con- 
fidence intervals of the parameters. The average of the 
squared semi-axes of the optimized ellipsoid, i.e. the ob- 
jective in ([T5]h is 0.07 times the average of the squared 
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FIG. 1. Projections and three dimensional ellipsoid with 
shadows of the linearized 95%-confidence regions before (light 
part) and after (dark part) the optimization. Computed with 
Pasveer parameters (|18|) . 



TABLE II. Radii of the confidence intervals before and after 
the optimization computed by (fT4|) with Pasveer parameters 
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FIG. 2. Comparison of current- volt age characteristics before 
and after the optimization for Pasveer parameters (|18|) . The 
continuous lines correspond to experiments with configuration 
C and T and the dashed lines correspond to experiments with 
optimal configuration £* and T* . 



(20) 



semi-axes of the ellipsoid not optimized. The results 
of the simulations before and after the optimization pro- 
cedure are shown in Fig. O The optimal configuration 
leads to higher current densities. A second example is 
taken from Mensfoort and Coehoorn, 2 who used 

p = (1.0 • 10~ 10 , 0.077, 4.25 • 10 26 ) , 

V = (0.1,...,0.9,l,2,...,10), 

C = (100,200,300), 

T = (235,270,305). 

Again we take boundaries like in ([T9]) and the resulting 
vectors are 

C = (50,187.8,296.8), 
T* = (200,274.7,350). 

Unlike the first example, the optimal lengths differ much 
more from the starting values and not all temperatures 
were raised, the lowest temperature was even reduced. In 
Fig. [3j we illustrate the confidence ellipsoid and the pro- 
jections and TABLE IIIII shows a comparison of the con- 
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FIG. 3. Projections and three dimensional ellipsoid with 
shadows of the linearized 95%-confidence regions before (light 
part) and after (dark part) the optimization. Computed with 
Coehoorn parameters ([20]) . 

fidence intervals in this case. The average of the squared 
semi- axes of the optimized ellipsoid here is 0.16 times the 
average of the squared semi-axes of the ellipsoid not op- 
timized, too. The results of the simulations before and 
after the optimization procedure are shown in Fig. SJ 
The optimized configuration leads to a wider range of 
current densities gaining more information. 



VI. CONCLUSION AND OUTLOOK 

We have applied optimum experimental design to or- 
ganic semiconductors modeled by the EGDM in order 
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TABLE III. Radii of the confidence intervals before and after 
the optimization computed by (|14|) with Coehoorn parameters 
(HQ)). 
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FIG. 4. Comparison of cur rent- volt age characteristics before 
and after the optimization for Pasveer parameters (|2Q|) . The 
continuous lines correspond to experiments with configuration 
C and T and the dashed lines correspond to experiments with 
optimal configuration £* and T*. 



to reduce the parameter uncertainty caused by mea- 
surement errors. The variances of the parameters, i.e. 
the average of the squared semi-axes of the confidence 
ellipsoid, were 0.07 and 0.16 times the average of the 
squared semi-axes of the not optimized confidence ellip- 
soid respectively. The proposed experiments, followed 
by a parameter estimation lead to a model, which is 
approximately 10 times more exact with respect to 
the assumed measurement error than the previous one. 
The classical methods for solving semiconductor models 
were extended to fit with the EGDM. The derivatives 
required for the optimization were computed via au- 
tomatic differentiation, leaving the error on machine 
precision level, even for higher derivatives. We used 
unipolar layer devices for simplicity, but the presented 
methods can also be applied to multi-layer devices, trap 
generation models and exciton rate equations or, more 



far-reaching, all other models which are based on the 
van Roosbroeck system. 
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